clear all;


spin_Delta=0.0;
spin_nvec=[0,0,1];
spin_nvec=spin_nvec/norm(spin_nvec);

sigx=[0,1;1,0];
sigy=[0,-i;i,0];
sigz=diag([1,-1]);
sig0=eye(2);

gen_SOC_term=@(bdvec) sqrt(2)*(i)*lambda_SOC*(bdvec(1)*sigx+bdvec(2)*sigy+bdvec(3)*sigz);


a1=[2.7508   -1.5882    4.4921];
a2=[0    3.1764    4.4921];
a3=[-2.7508   -1.5882    4.4921];

Rh_pos=[0.125000 0.125000 0.625000 ;
    0.625000 0.125000 0.125000 ;
    0.125000 0.625000 0.125000 ;
    0.125000 0.125000 0.125000];

Rh_pos=Rh_pos-0.125;

Rh_posC=Rh_pos;
Rh_posC(:,1)=Rh_pos(:,1)*a1(1)+Rh_pos(:,2)*a2(1)+Rh_pos(:,3)*a3(1);
Rh_posC(:,2)=Rh_pos(:,1)*a1(2)+Rh_pos(:,2)*a2(2)+Rh_pos(:,3)*a3(2);
Rh_posC(:,3)=Rh_pos(:,1)*a1(3)+Rh_pos(:,2)*a2(3)+Rh_pos(:,3)*a3(3);



Rh_posC



tetra1_posC=Rh_posC
tetra1_C=mean(tetra1_posC,1)

tetra2_posC=-Rh_posC
tetra2_C=mean(tetra2_posC,1)


newx_axis=zeros(4,3);
newy_axis=zeros(4,3);
newz_axis=zeros(4,3);
for ind=1:4
    zvec=tetra1_posC(ind,:)-tetra1_C;
    newz_axis(ind,:)=-zvec/norm(zvec);
end
newz_axis

newx_axis(4,:)=[1,0,0];
newx_axis(3,:)=[-1,0,0];
newx_axis(2,:)=[0.5,0.5*sqrt(3),0];
newx_axis(1,:)=[0.5,-0.5*sqrt(3),0];


for ind=1:4
    xvec=newx_axis(ind,:);
    zvec=newz_axis(ind,:);
    newy_axis(ind,:)=cross(zvec,xvec);
end




Rh_pos(2,:)+(Rh_pos(2,:)-Rh_pos(1,:))+[-1,0,1] %
Rh_pos(2,:)+(Rh_pos(2,:)-Rh_pos(3,:))+[-1,1,0]
Rh_pos(2,:)+(Rh_pos(2,:)-Rh_pos(4,:))+[-1,0,0]


Rh_pos(1,:)+(Rh_pos(1,:)-Rh_pos(2,:))+[1,0,-1]
Rh_pos(1,:)+(Rh_pos(1,:)-Rh_pos(3,:))+[0,1,-1]
Rh_pos(1,:)+(Rh_pos(1,:)-Rh_pos(4,:))+[0,0,-1]

Rh_pos(3,:)+(Rh_pos(3,:)-Rh_pos(1,:))+[0,-1,1]
Rh_pos(3,:)+(Rh_pos(3,:)-Rh_pos(2,:))+[1,-1,0]
Rh_pos(3,:)+(Rh_pos(3,:)-Rh_pos(4,:))+[0,-1,0]

Rh_pos(4,:)+(Rh_pos(4,:)-Rh_pos(1,:))+[0,0,1]
Rh_pos(4,:)+(Rh_pos(4,:)-Rh_pos(2,:))+[1,0,0]
Rh_pos(4,:)+(Rh_pos(4,:)-Rh_pos(3,:))+[0,1,0]


tetra1=Rh_posC;
tetra1C=mean(tetra1,1);

tetra2=Rh_posC;
tetra2(1,:)=tetra2(1,:)-a3;
tetra2(2,:)=tetra2(2,:)-a1;
tetra2(3,:)=tetra2(3,:)-a2;

tetra1
tetra1C=mean(tetra1,1)
tetra2
tetra2C=mean(tetra2,1)




figure(1);
clf;

linecc= [.5 .5 .5];
hold on;
% plot_tetra_v3(Rh_posC,0*a1,linecc,2)
% plot_tetra_v3(Rh_posC,a1,linecc,2)
% plot_tetra_v3(Rh_posC,a2,linecc,2)
% plot_tetra_v3(Rh_posC,a3,linecc,2)
% 
% %plot_tetra(-Rh_posC,0*a1,'r')
% plot_tetra_v3(-Rh_posC,a1,linecc,2)
 plot_tetra_v4(-Rh_posC,a2,linecc,2)
% plot_tetra_v3(-Rh_posC,a3,linecc,2)
% 
% plot_tetra_v3(-Rh_posC,a1+a2,linecc,2)
% plot_tetra_v3(-Rh_posC,a2+a3,linecc,2)
% plot_tetra_v3(-Rh_posC,a1+a3,linecc,2)



hexpts1=[];
hexpts1(1,:)=Rh_posC(3,:)+a1;
hexpts1(2,:)=Rh_posC(2,:)+a2;
hexpts1(3,:)=Rh_posC(1,:)+a2;

hexpts1(4,:)=Rh_posC(3,:)+a3;
hexpts1(5,:)=Rh_posC(2,:)+a3;
hexpts1(6,:)=Rh_posC(1,:)+a1;


hexpts2=[];
hexpts2(1,:)=Rh_posC(4,:)+a1;
hexpts2(2,:)=Rh_posC(1,:)+a1;
hexpts2(3,:)=Rh_posC(2,:)+a3;
hexpts2(4,:)=Rh_posC(4,:)+a3;
hexpts2(5,:)=Rh_posC(1,:);
hexpts2(6,:)=Rh_posC(2,:);



hexpts3=[];
hexpts3(1,:)=Rh_posC(4,:)+a2;
hexpts3(2,:)=Rh_posC(2,:)+a2;
hexpts3(3,:)=Rh_posC(3,:)+a1;
hexpts3(4,:)=Rh_posC(4,:)+a1;
hexpts3(5,:)=Rh_posC(2,:);
hexpts3(6,:)=Rh_posC(3,:);



hexpts4=[];
hexpts4(1,:)=Rh_posC(4,:)+a3;
hexpts4(2,:)=Rh_posC(3,:)+a3;
hexpts4(3,:)=Rh_posC(1,:)+a2;
hexpts4(4,:)=Rh_posC(4,:)+a2;
hexpts4(5,:)=Rh_posC(3,:);
hexpts4(6,:)=Rh_posC(1,:);



sig_x=[0,1;1,0];
sig_y=[0,-i;i,0];
sig_z=[1,0;0,-1];
sig_0=eye(2);

gen_spinh=@(nnvec) (sig_x*nnvec(1)+sig_y*nnvec(2)+sig_z*nnvec(3))/norm(nnvec);

eva_spinS =@(state_wf) [state_wf'*sig_x*state_wf,state_wf'*sig_y*state_wf,state_wf'*sig_z*state_wf];



newnvec_tetra=Rh_posC(3,:)-tetra1C;
newnvec_tetra=newnvec_tetra/norm(newnvec_tetra);


%%

vecsize=1.5;
arrowcolor='r';
arrowsize=55;
arrowwidth=10;



tripts=Rh_posC(1:3,1:3);
tripts_xaxis=newx_axis(1:3,1:3);

frame_x=newx_axis(1:3,1:3);
frame_y=newy_axis(1:3,1:3);


%quiver3(tripts(:,1),tripts(:,2),tripts(:,3),tripts_xaxis(:,1),tripts_xaxis(:,2),tripts_xaxis(:,3),'LineWidth',7,'Color','k')
for ind=3
    %plt_newD2_orbital(tripts(ind,:),tripts_xaxis(ind,:),frame_x(ind,:),frame_y(ind,:));
    quiver3(tripts(ind,1)-0.5*tripts_xaxis(ind,1)*vecsize,tripts(ind,2)-0.5*tripts_xaxis(ind,2)*vecsize,tripts(ind,3)-0.5*tripts_xaxis(ind,3)*vecsize,tripts_xaxis(ind,1)*vecsize,tripts_xaxis(ind,2)*vecsize,tripts_xaxis(ind,3)*vecsize,'LineWidth',arrowwidth,'Color',arrowcolor,'MaxHeadSize',arrowsize);

    %     quiver3(tripts(ind,1),tripts(ind,2),tripts(ind,3),frame_x(ind,1),frame_x(ind,2),frame_x(ind,3));
%     quiver3(tripts(ind,1),tripts(ind,2),tripts(ind,3),frame_y(ind,1),frame_y(ind,2),frame_y(ind,3));
% 
%     quiver3(tripts(ind,1),tripts(ind,2),tripts(ind,3),tripts_xaxis(ind,1),tripts_xaxis(ind,2),tripts_xaxis(ind,3),'LineWidth',7);
end





shiftvec=a2;
tripts=Rh_posC([1,2,4],1:3);
%tripts_xaxis=newx_axis([1,2,4],1:3);

newnvec=Rh_posC(3,:)-tetra1_C;
newnvec=newnvec/norm(newnvec);
newnvecx=[1,0,0];
newnvecy=cross(newnvec,newnvecx);
tripts_xaxis=zeros(3,3);
tripts_xaxis(1,:)=0.5*newnvecx-0.5*sqrt(3)*newnvecy;
tripts_xaxis(2,:)=0.5*newnvecx+0.5*sqrt(3)*newnvecy;
tripts_xaxis(3,:)=[-1,0,0];

tripts_xaxis=-tripts_xaxis;

frame_x=newx_axis([1,2,4],1:3);
frame_y=newy_axis([1,2,4],1:3);

for ind=3
    %plt_newD2_orbital(tripts(ind,:)+shiftvec,tripts_xaxis(ind,:),frame_x(ind,:),frame_y(ind,:));
    quiver3(tripts(ind,1)-0.5*tripts_xaxis(ind,1)*vecsize+shiftvec(1),tripts(ind,2)-0.5*tripts_xaxis(ind,2)*vecsize+shiftvec(2),tripts(ind,3)-0.5*tripts_xaxis(ind,3)*vecsize+shiftvec(3),tripts_xaxis(ind,1)*vecsize,tripts_xaxis(ind,2)*vecsize,tripts_xaxis(ind,3)*vecsize,'LineWidth',arrowwidth,'Color',arrowcolor,'MaxHeadSize',arrowsize);

 %     quiver3(tripts(ind,1)+shiftvec(1),tripts(ind,2)+shiftvec(2),tripts(ind,3)+shiftvec(3),frame_x(ind,1),frame_x(ind,2),frame_x(ind,3));
%     quiver3(tripts(ind,1)+shiftvec(1),tripts(ind,2)+shiftvec(2),tripts(ind,3)+shiftvec(3),frame_y(ind,1),frame_y(ind,2),frame_y(ind,3));
%     quiver3(tripts(ind,1)+shiftvec(1),tripts(ind,2)+shiftvec(2),tripts(ind,3)+shiftvec(3),tripts_xaxis(ind,1),tripts_xaxis(ind,2),tripts_xaxis(ind,3),'LineWidth',7);
end

for ind=1:3
    %[dot(tripts_xaxis(ind,:),frame_x(ind,:)),dot(tripts_xaxis(ind,:),frame_y(ind,:))]
end


%quiver3(tripts(:,1)+shiftvec(1),tripts(:,2)+shiftvec(2),tripts(:,3)+shiftvec(3),tripts_xaxis(:,1),tripts_xaxis(:,2),tripts_xaxis(:,3),'LineWidth',7,'Color','k')


R3mat=zeros(3);
R3mat(3,3)=1;
R3mat(1,1)=-0.5;
R3mat(2,2)=-0.5;
R3mat(1,2)=0.5*sqrt(3);
R3mat(2,1)=-0.5*sqrt(3);


tripts_xaxis=(R3mat*tripts_xaxis')';
shiftvec=a1;
tripts=Rh_posC([3,1,4],1:3);
%quiver3(tripts(:,1)+shiftvec(1),tripts(:,2)+shiftvec(2),tripts(:,3)+shiftvec(3),tripts_xaxis(:,1),tripts_xaxis(:,2),tripts_xaxis(:,3),'LineWidth',7,'Color','k')

frame_x=newx_axis([3,1,4],1:3);
frame_y=newy_axis([3,1,4],1:3);

for ind=1:3
    %plt_newD2_orbital(tripts(ind,:)+shiftvec,tripts_xaxis(ind,:),frame_x(ind,:),frame_y(ind,:));
    %quiver3(tripts(ind,1)-0.5*tripts_xaxis(ind,1)*vecsize+shiftvec(1),tripts(ind,2)-0.5*tripts_xaxis(ind,2)*vecsize+shiftvec(2),tripts(ind,3)-0.5*tripts_xaxis(ind,3)*vecsize+shiftvec(3),tripts_xaxis(ind,1)*vecsize,tripts_xaxis(ind,2)*vecsize,tripts_xaxis(ind,3)*vecsize,'LineWidth',arrowwidth,'Color',arrowcolor,'MaxHeadSize',arrowsize);

end


tripts_xaxis=(R3mat*tripts_xaxis')';
shiftvec=a3;
tripts=Rh_posC([2,3,4],1:3);

frame_x=newx_axis([2,3,4],1:3);
frame_y=newy_axis([2,3,4],1:3);
%quiver3(tripts(:,1)+shiftvec(1),tripts(:,2)+shiftvec(2),tripts(:,3)+shiftvec(3),tripts_xaxis(:,1),tripts_xaxis(:,2),tripts_xaxis(:,3),'LineWidth',7,'Color','k')
for ind=1:3
    %plt_newD2_orbital(tripts(ind,:)+shiftvec,tripts_xaxis(ind,:),frame_x(ind,:),frame_y(ind,:));
    %quiver3(tripts(ind,1)-0.5*tripts_xaxis(ind,1)*vecsize+shiftvec(1),tripts(ind,2)-0.5*tripts_xaxis(ind,2)*vecsize+shiftvec(2),tripts(ind,3)-0.5*tripts_xaxis(ind,3)*vecsize+shiftvec(3),tripts_xaxis(ind,1)*vecsize,tripts_xaxis(ind,2)*vecsize,tripts_xaxis(ind,3)*vecsize,'LineWidth',arrowwidth,'Color',arrowcolor,'MaxHeadSize',arrowsize);

end






% 
% 
% 
% %%
% tripts=Rh_posC(1:3,1:3);
% tripts_xaxis=newx_axis(1:3,1:3);
% 
% 
% vecsize=1.5;
% arrowcolor='r';
% arrowsize=55;
% arrowwidth=10;
% 
% %quiver3(tripts(:,1),tripts(:,2),tripts(:,3),tripts_xaxis(:,1),tripts_xaxis(:,2),tripts_xaxis(:,3),'LineWidth',7,'Color','k')
% for ind=1:3
%     quiver3(tripts(ind,1)-0.5*tripts_xaxis(ind,1)*vecsize,tripts(ind,2)-0.5*tripts_xaxis(ind,2)*vecsize,tripts(ind,3)-0.5*tripts_xaxis(ind,3)*vecsize,tripts_xaxis(ind,1)*vecsize,tripts_xaxis(ind,2)*vecsize,tripts_xaxis(ind,3)*vecsize,'LineWidth',arrowwidth,'Color',arrowcolor,'MaxHeadSize',arrowsize);
% 
% 
%     %quiver3(coorcenter(1),coorcenter(2),coorcenter(3),0,0,coordinatelen,'MaxHeadSize',15,'Color','k','LineWidth',7)
% 
%     %plt_porbital_convertD(tripts(ind,:),tripts_xaxis(ind,:),plt_xaxis(ind,:),plt_yaxis(ind,:));
% end
% 
% 
% 
% 
% shiftvec=a2;
% tripts=Rh_posC([1,2,4],1:3);
% tripts_xaxis=newx_axis([1,2,4],1:3);
% 
% newnvec=Rh_posC(3,:)-tetra1_C;
% newnvec=newnvec/norm(newnvec);
% newnvecx=[1,0,0];
% newnvecy=cross(newnvec,newnvecx);
% tripts_xaxis=zeros(3,3);
% tripts_xaxis(1,:)=0.5*newnvecx-0.5*sqrt(3)*newnvecy;
% tripts_xaxis(2,:)=0.5*newnvecx+0.5*sqrt(3)*newnvecy;
% tripts_xaxis(3,:)=[-1,0,0];
% 
% 
% 
% for ind=1:3
%     %plt_porbital_convertD(tripts(ind,:)+shiftvec,tripts_xaxis(ind,:),plt_xaxis(ind,:),plt_yaxis(ind,:));
% 
%     quiver3(tripts(ind,1)-0.5*tripts_xaxis(ind,1)*vecsize+shiftvec(1),tripts(ind,2)-0.5*tripts_xaxis(ind,2)*vecsize+shiftvec(2),tripts(ind,3)-0.5*tripts_xaxis(ind,3)*vecsize+shiftvec(3),tripts_xaxis(ind,1)*vecsize,tripts_xaxis(ind,2)*vecsize,tripts_xaxis(ind,3)*vecsize,'LineWidth',arrowwidth,'Color',arrowcolor,'MaxHeadSize',arrowsize);
% 
% 
% end
% 
% 
% 
% 
% 
% 
% %quiver3(tripts(:,1)+shiftvec(1),tripts(:,2)+shiftvec(2),tripts(:,3)+shiftvec(3),tripts_xaxis(:,1),tripts_xaxis(:,2),tripts_xaxis(:,3),'LineWidth',7,'Color','k')
% 
% % 
% % R3mat=zeros(3);
% % R3mat(3,3)=1;
% % R3mat(1,1)=-0.5;
% % R3mat(2,2)=-0.5;
% % R3mat(1,2)=0.5*sqrt(3);
% % R3mat(2,1)=-0.5*sqrt(3);
% 
% 
% %tripts_xaxis=newx_axis([3,1,4],:);
% 
% 
% R3mat=zeros(3);
% R3mat(3,3)=1;
% R3mat(1,1)=-0.5;
% R3mat(2,2)=-0.5;
% R3mat(1,2)=0.5*sqrt(3);
% R3mat(2,1)=-0.5*sqrt(3);
% 
% 
% tripts_xaxis=(R3mat*tripts_xaxis')';
% 
% 
% shiftvec=a1;
% tripts=Rh_posC([3,1,4],1:3);
% %quiver3(tripts(:,1)+shiftvec(1),tripts(:,2)+shiftvec(2),tripts(:,3)+shiftvec(3),tripts_xaxis(:,1),tripts_xaxis(:,2),tripts_xaxis(:,3),'LineWidth',7,'Color','k')
% 
% 
% for ind=1:3
%     %plt_porbital_convertD(tripts(ind,:)+shiftvec,tripts_xaxis(ind,:),plt_xaxis(ind,:),plt_yaxis(ind,:));
%     quiver3(tripts(ind,1)-0.5*tripts_xaxis(ind,1)*vecsize+shiftvec(1),tripts(ind,2)-0.5*tripts_xaxis(ind,2)*vecsize+shiftvec(2),tripts(ind,3)-0.5*tripts_xaxis(ind,3)*vecsize+shiftvec(3),tripts_xaxis(ind,1)*vecsize,tripts_xaxis(ind,2)*vecsize,tripts_xaxis(ind,3)*vecsize,'LineWidth',arrowwidth,'Color',arrowcolor,'MaxHeadSize',arrowsize);
% 
% 
% end
% 
% shiftvec=a3;
% tripts=Rh_posC([2,3,4],1:3);
% %tripts_xaxis=newx_axis([2,3,4],:);
% tripts_xaxis=(R3mat*tripts_xaxis')';
% %quiver3(tripts(:,1)+shiftvec(1),tripts(:,2)+shiftvec(2),tripts(:,3)+shiftvec(3),tripts_xaxis(:,1),tripts_xaxis(:,2),tripts_xaxis(:,3),'LineWidth',7,'Color','k')
% for ind=1:3
%     %plt_porbital_convertD(tripts(ind,:)+shiftvec,tripts_xaxis(ind,:),plt_xaxis(ind,:),plt_yaxis(ind,:));
%     quiver3(tripts(ind,1)-0.5*tripts_xaxis(ind,1)*vecsize+shiftvec(1),tripts(ind,2)-0.5*tripts_xaxis(ind,2)*vecsize+shiftvec(2),tripts(ind,3)-0.5*tripts_xaxis(ind,3)*vecsize+shiftvec(3),tripts_xaxis(ind,1)*vecsize,tripts_xaxis(ind,2)*vecsize,tripts_xaxis(ind,3)*vecsize,'LineWidth',arrowwidth,'Color',arrowcolor,'MaxHeadSize',arrowsize);
% 
% end






box on;
axis equal;






% plt_D0_orbital([0,0,1],[1,0,0],[0,1,1],-1);

%view([1,1,0.5])

%view([110,20])

view([20,30])



%camlight('headlight')


coorcenter=[-3.,-1.,2];
coordinatelen=1.5;
% quiver3(coorcenter(1),coorcenter(2),coorcenter(3),coordinatelen,0,0,'MaxHeadSize',15,'Color','k','LineWidth',7)
% quiver3(coorcenter(1),coorcenter(2),coorcenter(3),0,coordinatelen,0,'MaxHeadSize',15,'Color','k','LineWidth',7)
% quiver3(coorcenter(1),coorcenter(2),coorcenter(3),0,0,coordinatelen,'MaxHeadSize',15,'Color','k','LineWidth',7)



alltripts=[];

linecolor=[0.0,0.4,0.0];
linecolor='b'
linecolor=[0.0,0.4,0.0];
linecolor='cyan'
linewidth=7;

tri_pts=zeros(3,3);
tri_pts(1,:)=Rh_posC(1,:);
tri_pts(2,:)=Rh_posC(2,:);
tri_pts(3,:)=Rh_posC(3,:);

alltripts=[alltripts;tri_pts];
% pp=fill3(tri_pts(:,1),tri_pts(:,2),tri_pts(:,3),'b');pp.FaceAlpha=0.5;

%plot_polygon_line(tri_pts,linecolor,linewidth)

scatri_pts=tri_pts;
tri_ptsC=mean(tri_pts,1);
scatri_pts(:,1)=scatri_pts(:,1)-tri_ptsC(1);
scatri_pts(:,2)=scatri_pts(:,2)-tri_ptsC(2);
scatri_pts(:,3)=scatri_pts(:,3)-tri_ptsC(3);
scatri_pts=scatri_pts*0.25;
scatri_pts(:,1)=scatri_pts(:,1)+tri_ptsC(1);
scatri_pts(:,2)=scatri_pts(:,2)+tri_ptsC(2);
scatri_pts(:,3)=scatri_pts(:,3)+tri_ptsC(3);
%pp=fill3(scatri_pts(:,1),scatri_pts(:,2),scatri_pts(:,3),'k'); %pp.FaceAlpha=0.5;


tri_pts=zeros(3,3);
tri_pts(1,:)=Rh_posC(1,:)+a1;
tri_pts(2,:)=Rh_posC(3,:)+a1;
tri_pts(3,:)=Rh_posC(4,:)+a1;
alltripts=[alltripts;tri_pts];
% pp=fill3(tri_pts(:,1),tri_pts(:,2),tri_pts(:,3),'b');pp.FaceAlpha=0.5;
%plot_polygon_line(tri_pts,linecolor,linewidth)

scatri_pts=tri_pts;
tri_ptsC=mean(tri_pts,1);
scatri_pts(:,1)=scatri_pts(:,1)-tri_ptsC(1);
scatri_pts(:,2)=scatri_pts(:,2)-tri_ptsC(2);
scatri_pts(:,3)=scatri_pts(:,3)-tri_ptsC(3);
scatri_pts=scatri_pts*0.25;
scatri_pts(:,1)=scatri_pts(:,1)+tri_ptsC(1);
scatri_pts(:,2)=scatri_pts(:,2)+tri_ptsC(2);
scatri_pts(:,3)=scatri_pts(:,3)+tri_ptsC(3);
%pp=fill3(scatri_pts(:,1),scatri_pts(:,2),scatri_pts(:,3),'k'); %pp.FaceAlpha=0.5;



tri_pts=zeros(3,3);
tri_pts(1,:)=Rh_posC(1,:)+a2;
tri_pts(2,:)=Rh_posC(2,:)+a2;
tri_pts(3,:)=Rh_posC(4,:)+a2;
alltripts=[alltripts;tri_pts];
% pp=fill3(tri_pts(:,1),tri_pts(:,2),tri_pts(:,3),'b');pp.FaceAlpha=0.5;
%plot_polygon_line(tri_pts,linecolor,linewidth)


scatri_pts=tri_pts;
tri_ptsC=mean(tri_pts,1);
scatri_pts(:,1)=scatri_pts(:,1)-tri_ptsC(1);
scatri_pts(:,2)=scatri_pts(:,2)-tri_ptsC(2);
scatri_pts(:,3)=scatri_pts(:,3)-tri_ptsC(3);
scatri_pts=scatri_pts*0.25;
scatri_pts(:,1)=scatri_pts(:,1)+tri_ptsC(1);
scatri_pts(:,2)=scatri_pts(:,2)+tri_ptsC(2);
scatri_pts(:,3)=scatri_pts(:,3)+tri_ptsC(3);
%pp=fill3(scatri_pts(:,1),scatri_pts(:,2),scatri_pts(:,3),'k'); %pp.FaceAlpha=0.5;


tri_pts=zeros(3,3);
tri_pts(1,:)=Rh_posC(2,:)+a3;
tri_pts(2,:)=Rh_posC(3,:)+a3;
tri_pts(3,:)=Rh_posC(4,:)+a3;
alltripts=[alltripts;tri_pts];
%plot_polygon_line(tri_pts,linecolor,linewidth)

scatri_pts=tri_pts;
tri_ptsC=mean(tri_pts,1);
scatri_pts(:,1)=scatri_pts(:,1)-tri_ptsC(1);
scatri_pts(:,2)=scatri_pts(:,2)-tri_ptsC(2);
scatri_pts(:,3)=scatri_pts(:,3)-tri_ptsC(3);
scatri_pts=scatri_pts*0.25;
scatri_pts(:,1)=scatri_pts(:,1)+tri_ptsC(1);
scatri_pts(:,2)=scatri_pts(:,2)+tri_ptsC(2);
scatri_pts(:,3)=scatri_pts(:,3)+tri_ptsC(3);
%pp=fill3(scatri_pts(:,1),scatri_pts(:,2),scatri_pts(:,3),'k'); %pp.FaceAlpha=0.5;


linecolor=[0.0,0.8,0.0];
%linecolor='b'
linecolor='cyan';

pt1=Rh_posC(2,:);
pt2=Rh_posC(4,:)+a1;
%line([pt1(1),pt2(1)],[pt1(2),pt2(2)],[pt1(3),pt2(3)],'Color',linecolor,'LineWidth',linewidth)
pt1=Rh_posC(3,:);
pt2=Rh_posC(4,:)+a2;
line([pt1(1),pt2(1)],[pt1(2),pt2(2)],[pt1(3),pt2(3)],'Color',linecolor,'LineWidth',linewidth)

pt1=Rh_posC(1,:);
pt2=Rh_posC(4,:)+a3;
%line([pt1(1),pt2(1)],[pt1(2),pt2(2)],[pt1(3),pt2(3)],'Color',linecolor,'LineWidth',linewidth)

pt1=Rh_posC(3,:)+a1;
pt2=Rh_posC(2,:)+a2;
%line([pt1(1),pt2(1)],[pt1(2),pt2(2)],[pt1(3),pt2(3)],'Color',linecolor,'LineWidth',linewidth)

pt1=Rh_posC(1,:)+a2;
pt2=Rh_posC(3,:)+a3;
%line([pt1(1),pt2(1)],[pt1(2),pt2(2)],[pt1(3),pt2(3)],'Color',linecolor,'LineWidth',linewidth)


pt1=Rh_posC(2,:)+a3;
pt2=Rh_posC(1,:)+a1;
%line([pt1(1),pt2(1)],[pt1(2),pt2(2)],[pt1(3),pt2(3)],'Color',linecolor,'LineWidth',linewidth)




mirror_linewidth=5;
isor=1.5;
sphcenter=mean(alltripts,1)

sphere_scan=linspace(-2,2,51);
%# create coordinates
[xx,yy,zz] = meshgrid(sphere_scan,sphere_scan,sphere_scan+sphcenter(3));

%# calculate distance from center of the cube
rr = sqrt(xx.^2 + yy.^2 + (zz-sphcenter(3)).^2);
% p=patch(isosurface(xx,yy,zz,rr,isor));
% p.FaceColor = 'k';
% p.FaceAlpha = 0.15;
% p.EdgeColor = 'none';

ppval=0.15;

nvec=Rh_posC(2,:)-Rh_posC(1,:);
nvec=nvec/norm(nvec);

tstvec=[0,0,1];
cir_n1=cross(nvec,tstvec);
cir_n1=cir_n1/norm(cir_n1);
cir_n2=cross(cir_n1,nvec);

allphi=linspace(0,2*pi,121);
allphix=cos(allphi)*cir_n1(1)+sin(allphi)*cir_n2(1);
allphiy=cos(allphi)*cir_n1(2)+sin(allphi)*cir_n2(2);
allphiz=cos(allphi)*cir_n1(3)+sin(allphi)*cir_n2(3);

%plot3(isor*allphix(:),isor*allphiy(:),isor*allphiz(:)+sphcenter(3),'Color','k','LineWidth',mirror_linewidth)

%pp=fill3(isor*allphix(:),isor*allphiy(:),isor*allphiz(:)+sphcenter(3),'k'); pp.FaceAlpha=ppval;




nvec=Rh_posC(3,:)-Rh_posC(1,:);
nvec=nvec/norm(nvec);

tstvec=[0,0,1];
cir_n1=cross(nvec,tstvec);
cir_n1=cir_n1/norm(cir_n1);
cir_n2=cross(cir_n1,nvec);

allphi=linspace(0,2*pi,121);
allphix=cos(allphi)*cir_n1(1)+sin(allphi)*cir_n2(1);
allphiy=cos(allphi)*cir_n1(2)+sin(allphi)*cir_n2(2);
allphiz=cos(allphi)*cir_n1(3)+sin(allphi)*cir_n2(3);

%plot3(isor*allphix(:),isor*allphiy(:),isor*allphiz(:)+sphcenter(3),'Color','k','LineWidth',mirror_linewidth)
%pp=fill3(isor*allphix(:),isor*allphiy(:),isor*allphiz(:)+sphcenter(3),'k'); pp.FaceAlpha=ppval;






nvec=Rh_posC(4,:)-Rh_posC(1,:);
nvec=nvec/norm(nvec);

tstvec=[0,0,1];
cir_n1=cross(nvec,tstvec);
cir_n1=cir_n1/norm(cir_n1);
cir_n2=cross(cir_n1,nvec);

allphi=linspace(0,2*pi,121);
allphix=cos(allphi)*cir_n1(1)+sin(allphi)*cir_n2(1);
allphiy=cos(allphi)*cir_n1(2)+sin(allphi)*cir_n2(2);
allphiz=cos(allphi)*cir_n1(3)+sin(allphi)*cir_n2(3);

% plot3(isor*allphix(:),isor*allphiy(:),isor*allphiz(:)+sphcenter(3),'Color','k','LineWidth',mirror_linewidth)

%pp=fill3(isor*allphix(:),isor*allphiy(:),isor*allphiz(:)+sphcenter(3),'k'); pp.FaceAlpha=ppval;




nvec=Rh_posC(3,:)-Rh_posC(2,:);
nvec=nvec/norm(nvec);

tstvec=[0,0,1];
cir_n1=cross(nvec,tstvec);
cir_n1=cir_n1/norm(cir_n1);
cir_n2=cross(cir_n1,nvec);

allphi=linspace(0,2*pi,121);
allphix=cos(allphi)*cir_n1(1)+sin(allphi)*cir_n2(1);
allphiy=cos(allphi)*cir_n1(2)+sin(allphi)*cir_n2(2);
allphiz=cos(allphi)*cir_n1(3)+sin(allphi)*cir_n2(3);

% plot3(isor*allphix(:),isor*allphiy(:),isor*allphiz(:)+sphcenter(3),'Color','k','LineWidth',mirror_linewidth)
%pp=fill3(isor*allphix(:),isor*allphiy(:),isor*allphiz(:)+sphcenter(3),'k'); pp.FaceAlpha=ppval;




nvec=Rh_posC(4,:)-Rh_posC(2,:);
nvec=nvec/norm(nvec);

tstvec=[0,0,1];
cir_n1=cross(nvec,tstvec);
cir_n1=cir_n1/norm(cir_n1);
cir_n2=cross(cir_n1,nvec);

allphi=linspace(0,2*pi,121);
allphix=cos(allphi)*cir_n1(1)+sin(allphi)*cir_n2(1);
allphiy=cos(allphi)*cir_n1(2)+sin(allphi)*cir_n2(2);
allphiz=cos(allphi)*cir_n1(3)+sin(allphi)*cir_n2(3);

% plot3(isor*allphix(:),isor*allphiy(:),isor*allphiz(:)+sphcenter(3),'Color','k','LineWidth',mirror_linewidth)
%pp=fill3(isor*allphix(:),isor*allphiy(:),isor*allphiz(:)+sphcenter(3),'k'); pp.FaceAlpha=ppval;






nvec=Rh_posC(4,:)-Rh_posC(3,:);
nvec=nvec/norm(nvec);

tstvec=[0,0,1];
cir_n1=cross(nvec,tstvec);
cir_n1=cir_n1/norm(cir_n1);
cir_n2=cross(cir_n1,nvec);

allphi=linspace(0,2*pi,121);
allphix=cos(allphi)*cir_n1(1)+sin(allphi)*cir_n2(1);
allphiy=cos(allphi)*cir_n1(2)+sin(allphi)*cir_n2(2);
allphiz=cos(allphi)*cir_n1(3)+sin(allphi)*cir_n2(3);

% plot3(isor*allphix(:),isor*allphiy(:),isor*allphiz(:)+sphcenter(3),'Color','k','LineWidth',mirror_linewidth)
%pp=fill3(isor*allphix(:),isor*allphiy(:),isor*allphiz(:)+sphcenter(3),'k'); pp.FaceAlpha=ppval;


hold off;


box on;
axis equal;

view([0,30])
view([155,-35])

set(gca,'visible','off')
% camlight('headlight')
material dull

%%

[az,el] = view


%%

 %exportgraphics(gcf,'Pyrochlore_dxydx2y2_3dstate.eps','ContentType','vector','BackgroundColor','none')



